Initialisation

## [1] "number of cores available = 10"
#Phi[1] ; eta = valeur de fin
#Phi[2] = valeur du noeud
#Phi[3] = echelle
m <- function(t, eta, phi) (phi[1] + eta)/(1+exp((phi[2]-t)/phi[3]))
#=======================================#
param <- list(sigma2 = 0.05,
              rho2 = 0.1,
              mu = c(5,4,0.5),
              omega2 = c(0.5,0.1,0.01))

F. <- length(param$omega2) #dimension de phi
#=======================================#
t <- seq(2,6, length.out = 10) #value of times

# --- longitudinal data --- #
G = 40 ; ng = 4 #nombre de groupe et d'individu par groupe
n <- G*ng*length(t) ; N <- G*ng  #nombre total de data, et #nombre total d'individu

dt_NLME <- NLME_obs(G, ng, t, param, m)
Y <- dt_NLME$obs

SAEM avec simulation par MCMC

\[\log f(Y, Z ; \theta) = \langle \Phi(\theta) ; S(Y,Z) \rangle - \psi(\theta)\]

#Petite fonction pour retourner rapidement l'appel Phi[attr(Phi, i)] où i est '1', '2', ...,
`%a%` <- function(x,var){
  if(length(var)== 1) return(x[ attr(x,as.character(var)) ])
  lapply(var, function(v) x%a%v) %>% unlist
}

Phi <- function(sigma2, rho2, mu, omega2, nu2)
{
  p <-c(- n/(2*sigma2), -N/(2*rho2),
        - G/(2*omega2), G*mu/omega2 )
 
  return(p)
}

S <- function(eta, phi)
{
  s <- c(mean((Y - get_obs(m, dt_NLME, eta = eta, phi = phi) )^2 ), #1
         mean(eta^2),                                               #2
         apply(phi^2, 2, mean), apply(phi  , 2, mean) )             #3 et 4

  attr(s, '1') <- 1
  attr(s, '2') <- 2
  attr(s, '3') <- 2 + 1:ncol(phi)
  attr(s, '4') <- 2 +   ncol(phi) + 1:ncol(phi)
  return(s)
}

Metropolis Hastings

loglik.phi <- function(x, eta, Phi)
{
  s <- S(eta,x)
  sum( (Phi*s)%a%c(1,3,4) ) #indice de S_1 et S_{3,.} puis S_{4,.}
}
loglik.eta <- function(x, phi, Phi) 
{ 
  s <- S(x, phi)
  sum( (Phi*s)%a%c(1,2) )
}

SAEM

Initialisation

M <- 1 #nombre de simulation
u <- function(k) ifelse(k<200, 1, 1/(k-199) )

# ---  Initialisation des paramètres --- #
para <- param %>% sapply(function(x) x + runif(1))
para$rho2 = 0.2 ; para$mu <- c(6,3,1) ; para$omega2 <- rep(.1,3)

# --- Initialisation des chaines MC : Z_0 ---
Z <- list(eta = 1:M %>% lapply(function(i) rnorm(G*ng, 0, param$rho2)  %>% matrix(ncol = 1) ),
          phi = 1:M %>% lapply(function(i) matrix(rnorm(F.*G, param$mu, param$omega2), nrow = F.) %>% t ) )

Boucle

sim <- function(niter, Phih, eta, phi)
{
  M <- length(phi)

  eta <- 1:M %>% lapply( function(i)
    MH_High_Dim_para_future(niter, eta[[i]], sd = .036,                 loglik.eta, phi[[i]], Phih, cores = ncores-1 ))
  
  phi <- 1:M %>% lapply( function(i)
    MH_High_Dim_para_future(niter, phi[[i]], sd = c(.028, .036, .021), loglik.phi, eta[[i]], Phih, cores = ncores-1 ))
  
  list(eta = eta, phi = phi)
}

maxi <- function(S)
{
  list(sigma2 = S%a%1,
       rho2 =   S%a%2,
       mu =     S%a%4,
       omega2 = S%a%3 - (S%a%4)^2 )
}
sigma2 rho2 mu1 mu2 mu3 omega21 omega22 omega23
Oracle 0.0510249 0.094816 4.97881 3.950752 0.4798243 0.5195325 0.0975481 0.0072938
Initialisation 0.5404334 0.200000 6.00000 3.000000 1.0000000 0.1000000 0.1000000 0.1000000 
niter <- 50
MH.iter <- 10
res <- SAEM(niter, MH.iter, u, para, Phi, S, Z, sim, maxi, eps = 1e-3, verbatim = 2)
saveRDS(res, 'saem.rds')
gg <- plot(res)
## [1] "SAEM execution time = 15min 39sec"
Résultat de l’algo SAEM-MCMC
sigma2 rho2 mu.1 mu.2 mu.3 omega2.1 omega2.2 omega2.3
Valeur réelle 0.0500 0.1000 5.0000 4.0000 0.5000 0.5000 0.1000 0.0100
Valeur estimée 0.0525 0.2099 4.9736 3.9507 0.4825 0.2117 0.0945 0.0072
Rrmse 0.0491 1.0995 0.0053 0.0123 0.0350 0.5766 0.0548 0.2778